library(pacman)

### Load necessary packages
pacman::p_load( conflicted, data.table,tidyverse, ggplot2, dplyr, readr, tidyr, dplyr, Hmisc, stargazer, stringr, stringi, lfe, xtable, here, readxl,stringdist,mfx,sciplot) 

rm(list=ls())
conflict_prefer("first", "dplyr")
conflict_prefer("last", "dplyr")
conflict_prefer("here", "here")
conflict_prefer("list", "base")

specify_decimal <- function(x, k) format(as.numeric(round(x, k), nsmall=k))


load(here("..","..","Code","Replication","lev2019.Rda"))

###################################
####### MAIN TEXT          ########
###################################

#### Table 1

table(lev2019$treat)





#### Figure 1


lev2019_d<-lev2019[,list(outcome1=mean(outcome1,na.rm=TRUE),
                         outcome1_se=se(outcome1,na.rm=TRUE),
                         outcome3=mean(outcome3,na.rm=TRUE),
                         outcome3_se=se(outcome3,na.rm=TRUE)
                         
),by=c("treat")]
lev2019_d$treat<-as.character(lev2019_d$treat)
lev2019_d$treatname[lev2019_d$treat==0]<-"  Election\nCommissions\n\n(Control)"
lev2019_d$treatname[lev2019_d$treat==1]<-" Candidate\nFiltering\n\n(Treatment #1)"
lev2019_d$treatname[lev2019_d$treat==3]<-"Workplace\nMobilization\n\n(Treatment #3)"
lev2019_d$treatname[lev2019_d$treat==2]<-"Organizing\nCarousels\n\n(Treatment #2)"

ggplot(data=lev2019_d, aes(x=treatname, y=outcome1, fill=treatname)) + geom_bar(colour="black", stat="identity",position=position_dodge(),size=.3) + geom_errorbar(aes(ymin = outcome1-outcome1_se, ymax = outcome1+outcome1_se), width = 0.2)+ theme_bw()+ theme(legend.position="bottom",axis.text=element_text(size=16),axis.title=element_text(size=18),legend.text=element_text(size=12),plot.title = element_text(hjust = 0.5,size=22))+ coord_cartesian(ylim=c(2,5))+ geom_text(aes(label=specify_decimal(outcome1,2)), position=position_dodge(width=0.9), vjust=-1.3,size=5)+xlab("") + scale_fill_grey( name="",breaks=c("1", "2", "3"))+    xlab("") + ylab("Level of Anger\n")+ guides(fill=FALSE)+ggtitle("Panel A: Respondent Anger over Electoral Integrity")



ggplot(data=lev2019_d, aes(x=treatname, y=outcome3, fill=treatname)) + geom_bar(colour="black", stat="identity",position=position_dodge(),size=.3) + geom_errorbar(aes(ymin = outcome3-outcome3_se, ymax = outcome3+outcome3_se), width = 0.2) + theme_bw()+ theme(legend.position="bottom",axis.text=element_text(size=16),axis.title=element_text(size=18),legend.text=element_text(size=12),plot.title = element_text(hjust = 0.5,size=22))+ coord_cartesian(ylim=c(2,4))+ scale_y_continuous(breaks=c(2,3,4))+ geom_text(aes(label=specify_decimal(outcome3,2)), position=position_dodge(width=0.9), vjust=-1.3,size=5)+xlab("") + scale_fill_grey( name="",breaks=c("1", "2", "3"))+    xlab("") + ylab("Likelihood of Protesting\n")+ guides(fill=FALSE)+ggtitle("Panel B: Respondent Likelihood of Protesting")



#### Table 2

lev2019$treat_factor<-relevel(factor(lev2019$treat), ref = "0")

est1<-felm(outcome1~treat_factor|0|0|0, data=subset(lev2019))

est2<-felm(outcome1~treat_factor + male+logage+edu+econ+townsize+employed+voted2018|0|0|0, data=subset(lev2019))

lev2019$treat_factor<-relevel(factor(lev2019$treat), ref = "3")

est3<-felm(outcome1~treat_factor|0|0|0, data=subset(lev2019, treat!=0))

est4<-felm(outcome1~treat_factor + male+logage+edu+econ+townsize+employed+voted2018|0|0|0, data=subset(lev2019 , treat!=0))

lev2019$treat_factor<-relevel(factor(lev2019$treat), ref = "0")

est5<-felm(outcome3 ~ treat_factor, data=lev2019)

est6<-felm(outcome3 ~ treat_factor + male + logage + edu + econ + townsize + employed + voted2018, data=lev2019)

lev2019$treat_factor<-relevel(factor(lev2019$treat), ref = "3")

est7<-felm(outcome3~treat_factor|0|0|0, data=subset(lev2019, treat!=0))

est8<-felm(outcome3~treat_factor + male+logage+edu+econ+townsize+employed+voted2018|0|0|0, data=subset(lev2019 , treat!=0))


###################################
####### APPENDIX           ########
###################################


#### Table D1 

lev_summary<-lev2019[,list(male, logage,edu,econ,townsize,employed,voted2018)]

levtable<-stargazer(lev_summary,covariate.labels=c(
  "Male",
  "Age (log)",
  "Education",
  "Economic Situation",
  "City Size",
  "Employed",
  "Voted in 2018 Election"),header=FALSE,digits=3)


#### Tables D2 and D3

lev2019$treat_factor<-relevel(factor(lev2019$treat), ref = "0")

exp_sum<-lev2019[,list(male=mean(male,na.rm=TRUE),
                       logage=mean(logage,na.rm=TRUE),
                       edu=mean(edu,na.rm=TRUE),
                       econ=mean(econ,na.rm=TRUE),
                       townsize=mean(townsize,na.rm=TRUE),
                       employed=mean(employed,na.rm=TRUE),
                       voted2018=mean(voted2018,na.rm=TRUE) ),by=c("treat_factor")]
exp_sum$treat_factor<-as.character(exp_sum$treat_factor)
exp_sum<-exp_sum[order(treat_factor)]
exp_sum<-dcast(melt(exp_sum, id.vars = "treat_factor"), variable ~ treat_factor)
exp_sum$variable=NULL
rownames( exp_sum ) <- c("Male",
                         "Age (log)",
                         "Education",
                         "Economic Situation",
                         "Town Size",
                         "Employed",
                         "Voted in 2018 Election")
colnames( exp_sum ) <- c("Control",
                         "Treatment #1 (Filtering)",
                         "Treatment #2 (Carousels)",
                         "Treatment #3 (Workplace Mob.)")


#### Figures D1 and D2


lev2019_d<-lev2019[is.na(lev2019$putin_binary)==FALSE][,list(outcome1=mean(outcome1,na.rm=TRUE),
                                                             outcome1_se=se(outcome1,na.rm=TRUE),
                                                             outcome3=mean(outcome3,na.rm=TRUE),
                                                             outcome3_se=se(outcome3,na.rm=TRUE)
                                                             
                                                             
),by=c("treat","putin_binary")]

lev2019_d$treat<-as.character(lev2019_d$treat)
lev2019_d$treatname[lev2019_d$treat==0]<-"  Election\nCommissions\n\n(Control)"
lev2019_d$treatname[lev2019_d$treat==1]<-" Candidate\nFiltering\n\n(Treatment #1)"
lev2019_d$treatname[lev2019_d$treat==3]<-"Workplace\nMobilization\n\n(Treatment #3)"
lev2019_d$treatname[lev2019_d$treat==2]<-"Organizing\nCarousels\n\n(Treatment #2)"

lev2019_d_opp<-subset(lev2019_d, putin_binary==0)
lev2019_d_putin<-subset(lev2019_d, putin_binary==1)

ggplot(data=lev2019_d_opp, aes(x=treatname, y=outcome1, fill=treatname)) + geom_bar(colour="black", stat="identity",position=position_dodge(),size=.3)+ geom_errorbar(aes(ymin = outcome1-outcome1_se, ymax = outcome1+outcome1_se), width = 0.2) + theme_bw()+ theme(legend.position="bottom",axis.text=element_text(size=16),axis.title=element_text(size=18),legend.text=element_text(size=12),plot.title = element_text(hjust = 0.5,size=22))+ coord_cartesian(ylim=c(2,5))+ geom_text(aes(label=specify_decimal(outcome1,2)), position=position_dodge(width=0.9), vjust=-1.9,size=5)+xlab("") + scale_fill_grey( name="",breaks=c("1", "2", "3"))+    xlab("") + ylab("Level of Anger\n")+ guides(fill=FALSE)+ggtitle("Panel A: Respondent Anger over Electoral Integrity")

ggplot(data=lev2019_d_opp, aes(x=treatname, y=outcome3, fill=treatname)) + geom_bar(colour="black", stat="identity",position=position_dodge(),size=.3) + geom_errorbar(aes(ymin = outcome3-outcome3_se, ymax = outcome3+outcome3_se), width = 0.2) + theme_bw()+ theme(legend.position="bottom",axis.text=element_text(size=16),axis.title=element_text(size=18),legend.text=element_text(size=12),plot.title = element_text(hjust = 0.5,size=22))+ coord_cartesian(ylim=c(2,4))+ scale_y_continuous(breaks=c(2,3,4))+ geom_text(aes(label=specify_decimal(outcome3,2)), position=position_dodge(width=0.9), vjust=-1.9,size=5)+xlab("") + scale_fill_grey( name="",breaks=c("1", "2", "3"))+    xlab("") + ylab("Likelihood of Protesting\n")+ guides(fill=FALSE)+ggtitle("Panel B: Respondent Likelihood of Protesting")

ggplot(data=lev2019_d_putin, aes(x=treatname, y=outcome1, fill=treatname)) + geom_bar(colour="black", stat="identity",position=position_dodge(),size=.3)+ geom_errorbar(aes(ymin = outcome1-outcome1_se, ymax = outcome1+outcome1_se), width = 0.2) + theme_bw()+ theme(legend.position="bottom",axis.text=element_text(size=16),axis.title=element_text(size=18),legend.text=element_text(size=12),plot.title = element_text(hjust = 0.5,size=22))+ coord_cartesian(ylim=c(2,5))+ geom_text(aes(label=specify_decimal(outcome1,2)), position=position_dodge(width=0.9), vjust=-1.9,size=5)+xlab("") + scale_fill_grey( name="",breaks=c("1", "2", "3"))+    xlab("") + ylab("Level of Anger\n")+ guides(fill=FALSE)+ggtitle("Panel A: Respondent Anger over Electoral Integrity")

ggplot(data=lev2019_d_putin, aes(x=treatname, y=outcome3, fill=treatname)) + geom_bar(colour="black", stat="identity",position=position_dodge(),size=.3) + geom_errorbar(aes(ymin = outcome3-outcome3_se, ymax = outcome3+outcome3_se), width = 0.2) + theme_bw()+ theme(legend.position="bottom",axis.text=element_text(size=16),axis.title=element_text(size=18),legend.text=element_text(size=12),plot.title = element_text(hjust = 0.5,size=22))+ coord_cartesian(ylim=c(2,4))+ scale_y_continuous(breaks=c(2,3,4))+ geom_text(aes(label=specify_decimal(outcome3,2)), position=position_dodge(width=0.9), vjust=-1.9,size=5)+xlab("") + scale_fill_grey( name="",breaks=c("1", "2", "3"))+    xlab("") + ylab("Likelihood of Protesting\n")+ guides(fill=FALSE)+ggtitle("Panel B: Respondent Likelihood of Protesting")


#### Table D4 

est1<-felm(outcome3 ~ treat_factor, data=lev2019)

est2<-felm(outcome3 ~ treat_factor + male + logage + edu + econ + townsize + employed + voted2018, data=lev2019)

est3<-felm(outcome3 ~ outcome1 + treat_factor, data=lev2019)

est4<-felm(outcome3 ~ outcome1 + treat_factor + male + logage + edu + econ + townsize + employed + voted2018, data=lev2019)



#### Figure E3

load(here("..","..","Code","Replication","lev2018.Rda"))



lev2018_d<-lev2018[,list(list_outcome=mean(list_outcome,na.rm=TRUE),
                         list_outcome_se=se(list_outcome,na.rm=TRUE)
),by=c("treat")]
lev2018_d$treat<-as.character(lev2018_d$treat)
lev2018_d$treatname[lev2018_d$treat==0]<-"   No Activity\n\n\n(Pure Control)"
lev2018_d$treatname[lev2018_d$treat==1]<-"  Election\nCommissions\n\n(Control)"
lev2018_d$treatname[lev2018_d$treat==2]<-" Candidate\nFiltering\n\n(Treatment #1)"
lev2018_d$treatname[lev2018_d$treat==3]<-" Workplace\nMobilization\n\n(Treatment #2)"
lev2018_d$treatname[lev2018_d$treat==4]<-"Organizing\nCarousels\n\n(Treatment #3)"
lev2018_d$fill="4"
lev2018_d$fill[lev2018_d$treat==1]="3"
lev2018_d$fill[lev2018_d$treat==2]="2"
lev2018_d$fill[lev2018_d$treat>2]="1"

ggplot(data=lev2018_d, aes(x=treatname, y=list_outcome, fill=fill)) + geom_bar(colour="black", stat="identity",position=position_dodge(),size=.3) + geom_errorbar(aes(ymin = list_outcome-list_outcome_se, ymax = list_outcome+list_outcome_se), width = 0.2) + theme_bw()+ theme(legend.position="bottom",axis.text=element_text(size=14),axis.title.y=element_text(size=14),legend.text=element_text(size=12))+ coord_cartesian(ylim=c(1,5))+ geom_text(aes(label=specify_decimal(list_outcome,2)), position=position_dodge(width=0.9), vjust=-1.3,size=5)+xlab("") + ylab("Probability of Voting")+ scale_fill_grey( name="",breaks=c("1", "2", "3","4"))+    xlab("") + ylab("Likelihood of Voting")+ guides(fill=FALSE)


#### Table E7

lev2018$treat_factor<-relevel(factor(lev2018$treat), ref = "0")

est1<-felm(list_outcome~treat_factor |0|0|0, data=subset(lev2018))

est2<-felm(list_outcome~treat_factor + male+logage+edu+selfreportecon+citysize+employed+turnout|0|0|0, data=subset(lev2018))

lev2018$treat_factor<-relevel(factor(lev2018$treat), ref = "1")

est3<-felm(list_outcome~treat_factor|0|0|0, data=subset(lev2018, treat>0))

est4<-felm(list_outcome~treat_factor + male+logage+edu+selfreportecon+citysize+employed+turnout|0|0|0, data=subset(lev2018, treat>0))

lev2018$treat_factor<-relevel(factor(lev2018$treat), ref = "4")

est5<-felm(list_outcome~treat_factor |0|0|0, data=subset(lev2018, treat>1))

est6<-felm(list_outcome~treat_factor + male+logage+edu+selfreportecon+citysize+employed+turnout|0|0|0, data=subset(lev2018, treat>1))


#### Tables E8 and E9
lev2018$treat_factor<-relevel(factor(lev2018$treat), ref = "0")

exp_sum<-lev2018[,list(male=mean(male,na.rm=TRUE),
                       logage=mean(logage,na.rm=TRUE),
                       edu=mean(edu,na.rm=TRUE),
                       selfreportecon=mean(selfreportecon,na.rm=TRUE),
                       citysize=mean(citysize,na.rm=TRUE),
                       employed=mean(employed,na.rm=TRUE),
                       turnout=mean(turnout,na.rm=TRUE) ),by=c("treat_factor")]
exp_sum$treat_factor<-as.character(exp_sum$treat_factor)
exp_sum<-exp_sum[order(treat_factor)]
exp_sum<-dcast(melt(exp_sum, id.vars = "treat_factor"), variable ~ treat_factor)
exp_sum$variable=NULL
rownames( exp_sum ) <- c("Male",
                         "Age (log)",
                         "Education",
                         "Family Economic Situation",
                         "City Size",
                         "Employed",
                         "Voted in 2018 Election")
colnames( exp_sum ) <- c("Pure Control",
                         "Control",
                         "Treatment #1 (Filtering)",
                         "Treatment #2 (Workplace Mob.)",
                         "Treatment #3 (Ballot Box)")

